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ABSTRACT 

Incorporation of the effects of time-dependent diffusive propagation 
of galactic cosmic rays inside a modulating region whose basic parameters 
are slowly changing in time leads to a new prediction for the modulated 
density U(t) expected to be observed at a given time t, A first order 
perturbation analysis shows that if t:he expected density under 

completely stationary conditions at time t^, then the actual density under 
slowly varying conditions will be given by 

where t(K) is the average time spent by a particle of diffusion coefficient 
K inside the modulation region. An analysis of the behavior of t as a 
function of various modulating parameters in both an idealized one dimen- 
sional convective wind and a three dimensional radial wind shows that t 
can be greater than 100 days under reasonable values of these parameters, 

The general behavior is such that in a modulating region characterized by 
the distance to the modulating boundary R> the convective velocity V, and 
K, the average time t is proportional to R^/K in the limit of large K and 
R/V in the limit of small K for both geometries* This general behavior is not 



2 


appreciably affected by energy loss processes. Since t Is a function of K 
which is in turn a function of magnetic rigidity R and velocity g this 
model provides a natural physical explanation for observed rigidity depen- 
dent phase lags in modulated spectra sometimes referred to as cosmic ray 
"hysteresis". If all of the phase lag observed between 500 MeV protons 
and the Deep River Neutron intensity is attributed to the effects described 
here, the average distance to the modulating boundary during the last 
solar cycle is estimated to be - 55 a.u. 
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I. INTRODUCTION 

The dif f usion~convectioii model, first introduced by Parker (1958, 1963) 
to explain the 11 year modulation of galactic cosmic rays is based on the 
assumption of time stationary Interplanetary conditions. In the basic model, 
and essentially all variations of this model discussed to date, the solution 
for the modulated spectra at a particular epoch of the solar cycle is a 
function of the specific values assigned to the parameters characterizing 
the interplanetary medium Cl*®* solar wind velocity, diffusion coefficient, 
size of modulation cavity, etc) during that epoch (Joklpii, 1971, Fisk, 

1971) . The observed time variations in intensity and spectral shape over 
the whole solar cycle are then explained in terms of a gradual variation in 
these parameters. Such "quasi-steady" solutions have met with considerable 
success in explaining the observed time variations although some models have 
had to introduce parameters beyond those whose physical significance is im- 
plicit in the simple diffusion convection picture in order to fit all the 
observations as is discussed in some detail by Rygg, O' Gallagher, and Earl 
(1974). Only Parker (1965) has considered the time dependent propagation 
problem in modulation in any detail and he did not consider the effects on 
the expected modulated spectrum. Simpson (1964) considered the effect of 
changes in the modulating region which originate at the sun and are convected 
outwards but did not consider the diffusive propagation of the particles 
themselves. 

In this paper, a model is developed which incorporates to first order, 
the direct effects of the time dependent diffusive propagation of Inter- 
. stellar cosmic rays in a slowly changing Interplanetary medivim. Some con- 
cepts basic to this model but limited to a one dimensional convective re- 
gion were described in a preliminary report (O'Gallagher 1973). Here, 
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these concepts are developed more fully and extended to consideration of a 
three dimensional radial convective region. The model shows clearly that 
the effects of time dependent diffusive propagation can be quite significant. 
Furthermore the model predicts a rigidity dependent time delay or "lag" in 
the modulated spectra and as such may provide a natural, physically reason- 
able explanation for the so called "hysteresis effect". (Simpson, 1964; 
Balasubrahmanyan et al. 1968; Kane and Winckler 1969; 0 * Galliigher , 1969; 
Simpson and Wang, 1971, Rygg, et al. , 1974), A bonus of this model is that 
observed hysteresis effects, when interpreted in terms of the model, pro- 
vide a direct measure of the dynamical features of the modulating cavity 
which cannot be inferred from time-stationary models. The model is con- 
ceptually simple and it is not necessary to introduce any parameters beyond 
those implicit in the time-stationary model. 

The model is best Introduced in terms of the usual Fokker— Planck equa- 
tion for diffusive particle transport in the interplanetary medium. 




(oTu) 

ar9T 


( 1 ) 


Here U = U(r,t,T) is the particle density (of a particular species) with 
kinetic energy between T and T + dT, at heliocentric radius r and 
time t. 

K = K(3»Rft»t) is the effective interplanetary diffusion coefficient 
(which is here assumed to be a scalar, i.e. Isotropic diffusion) 
at r and t as a function of the particle velocity 8 and magnetic 
rigidity R. 

V is the solar wind velocity. 


a - 


is a factor which conq>ensates for the transition between 


T + 2To 
T + To 

relativistic and non-relativistic energy regions. 
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The terms involving derivatives with respect to T on the right hand side of 
equation 1 incorporate the so-called Compton-Getting effect (Gleeson and 
Axford, 1968 b) . These terms result from the transformation of the propagation’ 
equations from a frame moving with the radially diverging solar wind to a frame 
stationary in the solar system and in effect account for the effects of adiabatic 
energy losses in the expanding diffusing medium on the density spectrum observed 
in this stationary frame. 


The conventional treatment of equation 1 is to argue that in the long term 
modulation, changes in U with time are so small that 9U/9t 0 and the right hand 

side can be set identically equal to zero to obtain a solution. The original, 
solution of Parker {1963, 1965) has been modified to obtain approximate solutions 
including the effects of energy loss explicitly (Fisk and Axford 1969, Gleeson 
and Axford, I968) and sophisticated computer methods have been developed to 
obtain numerical solutions under a wide variety of assumptions in the inter- 
planetary medium (Fisk, 1971, Lezniak and Webber, 1971, Urch and Gleeson 1972). 

In all of these cases the solutions have been obtained under the assumption of 
stationary conditions. As we shall see, this assumption may not be strictly 
valid for some reasonable values of interplanetary medium parameters. 

II. OWE DIMENSIONAL MODEL 

To incorporate time dependent diffusive effects into the solution for 
modulated spectra even to first order it- is necessary in principle to solve 
equation (l) as it stands. Since the general solution of equation (l) cannot be 
obtained in closed form and even in special cases the solution is quite compli- 
cated it is most instructive to introduce the basic concepts in the context of a 
one dimensional limit. In this case the analog of equation (l) is 


— (x,t,T) = 


K 




9x 


( 2 ) 
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where x is the (one) spatial dimension along which particles propagate 
and K and V have been assumed independent of position for purposes of 
the model. 

It is important to realize in considering the time-dependent modu- 
lation problem, both here in the one dimensional case and later in three 
dimensions, that there are two distinct time-dependent processes involved; 
a) the time dependent diffusive propagation of the cosmic ray particles 
themselves and h) the time variations in the medium parameters (K and V) 
and boundary conditions. In the actual modulation dynamics, the effects of 
these processes are coupled inseparably. The conventional stationary solutions 
approximate a solution by ignoring the first process (a) entirely. The model 
introduced here provides a better approximation by l) holding the medivim 
parameters and boundary conditions constant and using equation (2) to deter- 
mine the time scale of the diffusive propagation process under a particular 
set of stationary conditions and 2) incorporating as a first order correction 
the effect of cosmic ray propagation with a non-zero diffusion time in a 
medium whose parameters are slowly changing. Implicitly, this approach 
assumes separability and neglects the coupling between processes a) and b) 
above. For instance, taking K and V independent of both time and position is 
clearly an approximation since, strictly speaking, any time variation in K • 
and/or V will propagate through the medium with velocity V producing a position 
variation at a given time. However since such changes are assumed to be slow 
in the formal analysis, the effects which would be produced by this coupling 
are small and of higher order. For instance, changes in the diffusive propa- 
gation time during propagation, or energy loss processes due to differential 
variation of the solar wind velocity are second order effects and neglected 
in this model. In effect then, K and V in equation (2) are parameters in a 
simple model in which the diffusion and convection processes are each repre- 
sented by a single quantity at a given time which is to be regarded as a 
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characteristic or average value throughout the modulation region. Although 
not usually discussed in detail, this interpretation of K and V is exactly 
the interpretation given to modulation parameters in virtually all of the 
usual stationary treatments. 

In the one dimensional case, since the volume element of the diffusive 
medium is not expanding, individual particle energies remain unaffected and 
the energy charge terms do not appear. Thus equation (2) is similar to the 
"classical" diffusion-convection modulation equation neglecting adiabatic 
energy loss effects (Parker, 1963). This one dimensional analog of solar 
modulation is represented schematically in Figure 1. A diffusive medium is 
convected past an observer at x = 0 with velocity V to a distance X where 
there is a boundary beyond which particles move freely without scattering. 

For purposes of the model the solution to equation (2) is well represented 
in terms of the distribution function P(t), which is the probability per unit 
time for finding a particle (A) which crosses the boundary (x = X) at t = 0, 
inside the bo\mdary (at x<X) between t(>0) and t+dt- Parker (1965) has 
analyzed this problem in some detail and his solution for P at the observer 
(x = 0) can be expressed with a slight change in notation as 


y - (Vt+X)^ 

P(x,V,K,t) = ^ UKt 

(WKt)\ 


(3) • 


Equation (3) is essentially the Green's function solution for equation (2), 
In analogy to the real physical modulation problem, the particle density at 
and beyond the boTondary is taken to be a constant, U^, for all time. Then 
the density observed at x = 0 at a given time t^ is 

u(o,t^) = 


£ 


P[X,(t^-t),V,K] dt 
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Strictly speaking, equation (U) is valid only in the limit that X, V and K 
are all constant with respect to time in which case integration of the right 
hand side of equation (U) (see Appendix B) yields the familiar stationary 
solution 

_VX 

U (X,V,K) = U e ^ (5) 

s’’ o 

If in fact K, X and V.are slowly varying in time, then Equation (U) together 
with equation (3) provides the basis for incorporating the effects of these 
variations. The effect of variations in K, X and V are coupled to the rate 
at which changes are occuring, so a perturbation approach has been carried 
out. The detailed procedure followed is described in Appendix A and the 
result is discussed in some detail below. However at this point it is 
appropriate to emphasize the physical basis for the effect produced. This 
can be seen very simply from the distribution of diffusive propagation times 
from equation 3 plotted in Figure 2 for U different sets of interplanetary 
(one-dimensional) conditions. The distributions are plotted as P(t^^t), a 
function of the time in the "past" at which a particle seen now (at 't=t^) 
crossed the boundary. As a means of characterizing the time scales of the 
distributions, the time delay, t^, for which the distribution is a maximum 
is indicated for each case. The fundamental point is that the characteristic 
time scales depends on K (as well as X and V) which is a function of R and 
3. This implies that an observed spectrum of particles having a range of R 
and 3 will have sampled the characteristics of the interplanetary medium over 
a range of time scales in the past. As these characteristics change slowly 
this produces a rigidity dependent delay in the modulation at some rigidities 
with respect to other rigidities. 
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A particular exeimple of the dependence of the time scale on K is 

shown by curves (a) and (b) in Figure 2, which have the same values of 

V(V^ = - 300 km/sec) and X(X^ = = 1+0 a,u. ) but different K(K^ - 2K^ 

= 4 X 10^^ cm^/sec). Note first that the most probable time t^ is greater 

than, 100 days in both cases for values of X and V which are not untypical 

of those in conventional models and a value of K appropriate to ^-l<-00 MeV 

protons based on observed magnetic field power spectra (Jokipii and Coleman 

1968 ). Furthermore note that a difference of a factor of 2 in K corresponds 

to a difference of more than Uo days in the most probable delay time during 

which the characteristics of the medium may change slightly and therefore 

require a correction to the predicted spectrum. In the formal analysis , the 

first order effect of such corrections will depend on the average time x 

which particle spends in the diffusive medium, rather than the most probably 

time t , but the basic behavior of the time scale is well Illustrated by t . 
m m 

The parameters X and V of course do not depend on particle parameters 
and so do not directly affect the observed spectrum in the same way that K 
does. However they do affect the relevant time scales in a very direct way. 
Differences in X are illustrated by curves b) and d) which have the same V 
(V^ = = 300 km/sec) and K(K^ = = 4 x 10^1 cm^/sec) but have different 

X{X^ = Xj^/2 = 20 a.u. ). Particles from the boundary reach the observer much 
faster in a smaller region so that the characteristic time scale is strongly 
dependent on X. Differences in V are illustrated by curves (c) and (d) which 
have the same X(X^ = X^ = 20 a.u. ) and K(K^ = — 4 x 10^^ cm^/sec) but have 

different V(V^ = 2V^ = 6 OO km/sec). The larger value of V in case c) causes 
the exponential decrease of P(t) to set in earlier and shifts 

the most probable time to shorter values of (t^— t). The quantitative depen- 
dence of the time scales on all three parameters will be evaluated in §IV. 
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All of the curves in Figure 2 are normalized so that the integral of 
P(t^-t) for all t < t^ is equal to the protability that the particle actually 
arrives at x = 0 without being convected back out beyond x = X, or simply 

^vx/k 

e . Note in particular that curves b and c have the same value of VX/K 

although the differences in the individual parameters give rise to rather 
different time scales. 

Although the function P(X,V,K,t) given by equation (3) and plotted in 
Figure 2 serves to illustrate the dependence of the time scales on X,V, and 
K) it is not a completely valid one dimensional physical analog of modula- 
tion in that the diffusive-convective medium is assumed to be infinite in 
extent on the -x direction and in certain limits this leads to non-physical 
results so that the diffusing medium must be bounded. For simplicity in this 
model, the effects of diffusion-convection in a finite region, assumed to be 
roughly symmetric, can be well approximated formally by assiuning the distri- 
bution is characterized by a scale length equal to X. If the origin is 
redefined as the center of such a symmetric region then the resultant behavior 
at the origin is similar to that given by equation 3 but modified by the exponential 
decay factor e ' resulting from the gradual escape from the finite region 
yeilding the approximate expression for P(t) 


P(X,X«X,V,K,t) ^ 2 ^-3/2 

{UTfK)l/2 


exp 


(vt + x)^ 

UKt 



( 6 ) 


The effect of such a gradual escape is illustrated in Figure 2 by the dashed 
lines where we now have the additional requirement that the position of the 
observer at X must be such that JC«X. The exact solution of equationC2), for a 
finite region with free escape boundaries at x = ix is derived in an Appendix 
(c) and yield.s, a more complex expression for P(Xpc ,V,K,t) given by 
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V(X-x) „ 

P(X,x,V,K,t) = — e cos e 

^ n=0 

■faut the "basic physical concepts and characteristic time scales corresponding to 
equation (j) are essentially identical to those for equation (6) as will be shown. 

for the perturbation analysis it is most convenient to characterize the 
propagation delay by the average propagation time defined as follows: 


(n+^)^TT^Kt ^ V^tl 

X2 4K J 

(T) 


^ ^ ^ tP(X.V.K.t)dt 

^ ^“p(X,V,V,Kt)dt 

In terms of this parameter it is shown in Appendix A that the solution for the 
modulated density U(t^) corrected to first order for variations in modulating 
parameters is 


U(t ) = U (t ) - T 
o s o 


dt 



( 9 ) 


were U^(t^) is the solution obtained under stationary conditions using modulating 
parameters evaluated at t = Furthermore if dU^/dt is roughly constant over a 

time T the results of the analysis can be written in the simple and physically 
sensible form 


U(t^) [t^ - t(k)] (10) 

where we have explicitly noted that t is a function. of K and therefore R and 
&, At large rigidities, K becomes large and t(K) is small so that the stationary 

solution is quite accurate. At lower rigidities x can range from a few days to 
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several months with reasonable parameters in this one dimensional model so 
that equation (lO) provides a completely natural qualitative explanation for 
what is commonly referred to as the "hysteresis effect". 

The correction for these time varying effects is thus reduced to a 
determination of t(K) which can be accomplished relatively simply for one 
dimensional propagation using equations (6) or (?) and (8). Discussion of 
the detailed results of this determination are deferred to §IV where the I 

three dimensional problem can be discussed at the same time. 

Ill, A THREE DIMEarSIOMAL RADIAL WIMP 

Although the mathematical formulation involved for modulation in a 
three dimensional solar cavity is considerably more complicated, the basic 
physical concepts and qualitative solutions are identical to those intro- 
duced in the one dimensional treatment. If one assumes isotropic diffusion 
in a sperically symmetric, hcanogeneous interplanetary medium, equation (l) 
takes the form 


3U ^ K 9 / 2 

3t = 17 > 


v_ 

r2 


17 * 


2V ±_ 
3r 9T 


(aTU) 


( 11 ) 


A similar problem has been solved analytically for a special case of anisotropic 
diffusion subject to boiuidary and initial conditions appropriate for solar flare 
particles (Lupton and Stone, 1973) but equation (ll) has not been solved 
analytically for the general case with boundary conditions appropriate for 
solar modulation. 

There are two aspects in which equation (ll) differs from equation (2) 
and which make its solution more difficult. These are 

a) The three dimensional diffusive propagation described by the 
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first two terras on the right hand side 

h) The adiabatic energy losses described by the last term on the right. 
Since the diffusive process depends on K which in turn depends on energy, 
these two processes are coupled in a way that precludes analytic solution. 
Therefore it is necessary to analyze the spatial propagation process separa- 
tely from the energy loss process. This is accomplished by assuming that 
the last term in equation (ll) does not appreciably affect K in the analyses 
of the spatial transport problem. This of course is not true at all energies 
and the effects of the breakdown of this assumption on the results obtained 
will be considered in Section IV. 


Using notation similar to that in Section II the Green's function solution 
to equation (ll) for the spatial probability density in a radial wind for 
finding a particle which crosses a modulating boundary at r = R at t = 0 
and diffuses with constant K inside the boundary at r < R and t > 0 has been 


shown in a slightly modified notation by Parker {I965) to be of the form. 


P (R,r,K,V,t) = Z 
n=l 


a Q(w , 
n n 


-w V^t 
n 

K 


where Q{w^,Vr/K) are particular solutions chosen such that P(R,R,V,K,t>0) = 0 
and the a^ are constants chosen to satisfy the initial conditions. This is 
analogous to the one dimensional solution for P(x,x,V,K,t) given by equation (7). 
Assuiming that the series for P^ converges the conclusions are the same as for 
the one dimensional case. That is, one defines an average propagation time as 


T 

r 


^ tp(R,Y,K,r,t) dt 
P(R,V,K,r,t) dt 


( 13 ) 



Ik 


and the perturbation analysis in Appendix A leads to 


U(r,to) = - T — [^^(r.t^)] 


't = to 


dU 

or if — = constant, 
dT 


U(r.t^) ^ Ujr,t^ - T^] 


ilk) 


(15) 


The requirement that dU /dt he roughly constant means that equation (l5) will 

s 

not he valid for a time interval x near the reversal in phase of the modu- 
lation when higher order effects will need to he included. But for the long 
periods between solar maximum and solar minimum this is a general result 
which applies to all stationary models including the most sophisticated 
computer derived solutions which include the effects of energy loss. As with 
the one dimensional case the problem reduces to finding x hut here the analysis 
is complicated hoth hy the radial geometry and energy loss effects since K is 
not strictly constant through the propagation. Clearly however, 

under normal assumptions of spherical synometry and modulation 
in a radial wind one must expect non-negligible diffusive propagation times 
which depend on K and therefore R and 3. If these times are sufficiently 
long they will produce "phase lag" effect on the observed spectrum given by 
equation ( 15 ), which is qualitatively similar to the observed so called 
(hysteresis effect". In the following sections, the evaluation of x(K) in 
both the one and three dimensional models is discussed in some detail. 
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IV. EVALUATION OF AVERAGE DIFFUSIVE PROPAGATION TIMES 

The behavior of the average propagation time in one dimension can 
be determined in a straight forward manner for a given P(t) using equation 
(8). Consider first the behavior including the effects of escape from a 
finite region as approximated by equation (6). After substitution in equa- 
tion (8) both the numerator and denominator can be integrated as described 
in Appendix B yielding 

t^(X,V,K) = 

This function is plotted in Figure 3 to illustrate the predicted dependence 
of on K for a fixed convective velocity V = 300 km/sec and a range of 
values of X. Next consider the average time delay based on the exact solution 
given by equation (7) as described in Appendix C for one dimensional diffusion 
convection in a symmetric boimded region for an observer near the origin. 
Substitution of P(t) from equation (7) into equation (8) and integrating 
yields 

I (n+i) (-1)*^ [(n+i)^ 7 t2 + 

T (x,x«x,v,K) = ^£2 (i"?) 

^ S (n+|) (-1)^ [(n+i)2 u2 + (||)V 
n=0 

Both the nimerator and denominator of equation (17) are alternating series 
which converge. Values of calculated for the first 50 terms in each 
series are indicated in Figure 3 by x*s. The agreement between the exact 
solution and that based on approximate representation is extremely good. 
Therefore for all intents and purposes the analytical form given by 
equation (l6) can be used to calculate x for any desired values of X,V, 


x2/K 

V2/K + UK/X^ 




(l6) 



l6 


and K rather than the much more cumiberson representation given by equa- 
tion (it)* 

There are two limiting cases of special interest; Case a) the diffusion 
limit (V O) and Case b> the convection limit (K O). Both equations (l6) 
and ( 17 ) give identical results in these limits. In particular if we define 
as the propagation time in the diffusion limit and as the time in the 
convection limit we find 

Tp(X,K,V^) = g ( 18 ) 

and 

t^(X,K-K),V) =1 . (19) 

This latter result, which is surprisingly simple is a consequence of the fact 
already mentioned that as V becomes large compared to 2K/X the exponential 
decrease with time for P(X,V,K,t) in equation (3) prevents there being any 
substantial contribution to the integral in equation (8) from times much 
longer than X/V. Looked at more physically the result expressed in equation 
(19) is a manifestation of the fact that although smaller diffusion coeffi- 
cients will decrease the diffusive propagation speed (and this might be 
thought to increase the time) they also decrease the probability that a given 
particle will reach the observer at x = 0 before being convected back across 

— vx /k 

the boundary* Since we are interested only in those fraction (e~ ) of the 

original particles which are observed at x = 0, these processes set an upper 
limit on the average time that such an observed particle can spend inside 
the boundary. This upper limit as indicated by equation ( 19 ) is simply 
equivalent to the time required to convect a given element of the diffusive 
medium through the distance from the observer to the boundary of the medium. 
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Note finsdly that equation (l6) which provides an exceedingly good approxi 
mation to the exact solution can he written simply as 



The average propagation time in a 3 dimensional, radial convective 
medium cannot he evaluated for all K and V as was possible in one dimension. 
Therefore for purposes of this paper, to illustrate the physical validity 
of the model and to provide a quantitative approximation to the expected 
behavior we have considered the two limiting cases discussed in the one 
dimensional propagation mode above. In particular in the diffusion limit 
of Case a) it is shown in Appendix C that is given by 

(R,r,V-^0,K) = (21) 

r 

In deriving equation (2l) the additional condition that r << R has been imposed 
so that this result will be a good approximation at the orbit of earth only for 
R > 5-10 a.u. In the convection limit (Case b) , it is shown in, Appendix C 
that the propagation equation (ll) approaches the form of the one dimensional 
propagation equation in the limit K « Vr. Therefore from equation (l6) or 
(17) we find 

Tjj(R.r,V,K-vO) = ^ (22) 

- V 

which is again an upper limit equal to the characteristic convection time. 

Not only in the stated limit is equation (22) a very accurate approximation 
to the three dimensional behavior but since t is based on an average over 
all r < R, it will be a good approximation as long as K « VH. Thus equa- 
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tions (21) and (22) represent valid approximate solutions for in the limits 
of cases (a) and (b) analogous to the cases for the one dimensional mode above. 

To approximate the behavior of in the transition region between the 
two limiting cases we have used a form identical to that given by equation 
(20) which has been shown in the one dimensional case to be virtually 
identical to that given by the exact solution. The approximate behavior of 
t^(V,R,K) is plotted in Figure It as 


t^(R,V,K) 



" X 

H 

+ 

Hj 

1 L, 

-1/2 

1 r 2 /k 


[■^Dr" ^Cr^ 


V2/K + 36K/R2 


1/2 


(23) 


which gives precisely the correct behavior in the limiting regions (solid 
lines) and smoothly connects these limits through a transistion region (dashed 
lines) where Vr < K < VR where it should be a good approximation. From 
equation (23) it is apparent that the characteristic break, between the two 

VR 

limiting regimes occurs at K 'v ^ . The above discussion shows that the 
physical concepts illustrated in the one dimensional model are directly 
applicable to propagation in a radial, wind with only small changes in the 
quantitative behavior to be expected. Again we have an R^/K dependence of 
T at large K. As K becomes smaller this increases until it reaches the limit 
set by the convection time. Delays greater than 100 days are seen to be 
predicted for values of R, V and K which are not unreasonable and in fact are 
commonly used in other stationary solutions of the modulation problem. Clearly 
the effect of these delays cannot be neglected at low rigidities (small K) and 
Figure 4 and equation (15) provides a direct way to incorporate these correc- 
tions into any stationary model. 

Finally we consider again the effect of energy loss processes from the 
last term in equation (ll). Effectively what we have done in obtaining the 
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results in Figure k is to determine the approximate propagation time as 

if a given particle diffused in from r = R while maintaining the same value 

of K throughout the process. Clearly if the particle is losing energy 

during the process and K is in some way dependent on energy this assumption 

cannot be completely accurate. But it is easily shown that the inclusion 

of such effects cannot have any appreciable effect on the behavior of t (K) 

r 

as approximated by the curves in Figure U. Consider first that the energy 

loss process is also characterized by a time scale t. given by Parker 

J.OS s 

(1965) as 

2 a(T)V (2J4) 

~ K 

" V 

to a good approximation. This time is defined such that the energy E(t) of 

a particle which crosses the boTjndary at t = 0 with E( 0 ) = E will be 

o 

E(t) = E^ exp ■ (25) 

Note that = t^(R,V,K « VR/ 2 ) which is the upper limit on the propagation 

time. Thus, effectively by definition, the delay time in the large K limit 
will be a fraction of the characteristic energy loss time so that energy loss 
effects will be small and the basic assumption of constant K remains valid. 

At low energies since K is a decreasing function of kinetic energy (Jokipii 
and Coleman, I965) down to at least a few tens of MeV the propagation delay will 
increase and approach the characteristic loss time In this limit, energy 

losses are increasingly important but the propagation delay is constant 
independent of K. Thus it will remain unaffected by the energy loss processes 
and the approximate behavior in Figure U will not be modified appreciably. 

This can be seen more formally by representing the average effect of 
the energy loss process in terms of the propagation of a particle whose 
diffusion co-efficient is a monotonlcally decreasing function of time during 
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diffusion. Specifically consider the behavior of a particle which crosses 

the boundary at r = R at t = 0 with K = and is observed at r<B at t>0 

with K = K after undergoing propagation with a diffusion co-efficient K (t) 
r ^ 

such that K >K (t)<K at all times. First compare the average time t for 
Or— p -- r ^ 

two idealized cases: 

a) Propagation with constant diffusion coefficient so that 

T is given "by equations (21-23) and 
r 

b) propagation with Kp(t) = for t<t^ and Kp(t) = 

t>t^. The behavior in the two cases is identical for t<t^. In particular 
at t = t^ the expected intermediate distribution in position r' P(R,r",K^,V,tj^) 
is identical. Furthermore a lower limit for the distribution in propagation times 
from any r' to the observer at r is given by p(r^,r ,Kp,v,t -t^) . Thus a lower limit to 

the average time from any r'' to r is given by equations (21-23) (with substitution 

of r"for R and t-t-j^ for t). Substitution of Into equations (21-23) shows 

that the average propagation time from r' to r is greater in Case (b) than in 

Case (a) for all r". From this it follows that t(Kj^) > where 

for any time during the diffusion process. Thus for monotorlcally decreasing 

K (t) with time we 
P 

. conclude that in general x[Kp(t) < K^] > Conversely 

if K (t) > K it follows that tIK (t) > K ] i t(K ). If we apply these 
p r p — r X 

inequalities to particles which enter a modulating region with a spectrum 

in K and are observed with a spectrum in K we have the following cases 
o ^ 

of interest: 

1) K In this case we apply equation (2l) for K and find that 

' r D , r 

T <<T • Thus we know that K ^ K - K = constant so that t K^/6 k 

r loss r o p r 

is accurate. 

2) K «VR/6. In this case K <K we apply equation (22) and find that 

0 r o 

t(K ) ~ t(K ) and by the inequalities above we must have T(Kp) 

o V r V 
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Thus T = t(K ) as before is accurate in this limit also. 

T 

3) VR/6 and VR/6 (with of course K^<K^). Application of 

equation (20) and the inequalities above in this transition region yields 
^ ^ "^(Kp) ^ so that T = t(K^) is a reasonable approximation. 

In all three cases the propagation delays obtained for a given K 

r 

without considering the effects of energy loss are not appreciably altered 
by inclusion of these effects. It must be emphasized however that an 
accurate application of the entire model should include energy loss effects 
in calculating the stationary solutions in equation (l3). 

V. SUMMARY AMD DISCUSSION 

In the preceding sections of this paper, the effect of time dependent 
diffusive propagation in a medium whose characteristics are slowly changing, 
has been examined for the first time. These effects are found to be appre- 
ciable for some reasonable values of modulation parameters. A model which 
incorporates these effects as a first order perturbation on the usual time 
independent solution has been developed. The model predicts a "phase lag" 
between the response of cosmic ray intensities at different rigidities to 
changes in modulation parameters which is strikingly similar to observed 
hysteresis effects. In brief the most important conclusions of this analysis 
as they apply to modulation in a radial convective region are as follows: 

(1) The effects of time dependent diffusion propagation can be formally 
analyzed in terms of the average time spent by a particle inside 
the modulation region. 

(2) This average time is a function of the particle diffusion coefficient K. 

(3) The solution for modulated density of cosmic rays observed at time t 
corrected to first order for variations in modulation parameters is 

U(t) - U^(t - T (K)) 

& i 
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where U is the density predicted by the stationary solution for K. 
s 

(k) In the limit of large diffusion coefficients x^(K) is inversely pro- 
portional. to K and directly proportional to 

(5) In the limit of small diffusion coefficients the propagation time, 
becomes ec^ual to R/V, which is an upper limit on independent of K 

n 

as long as K « VR/6. 

(6) Propagation times the order of 100 days or more are consistent with 
reasonable values of V,R, and K. 

(T) This behavior is not appreciably affected by the action of adiabatic 
energy loss processes. 

In recent years much effort has been spent in developing modulation 
models which can explain the so called "hysteresis effects". Considerable 
success has been achieved by such authors as Van Hollebeke et al. (1972, 
1973), Burger and Swannenburg (l973) Bedijn et al (1973) by introducing 
changes in the rigidity dependence of K or incorporating rigidity dependence 
into other time variables parameters such as R. Processes which could be 
described by such phenomenological models may in fact be operative but 
there is as yet no independent evidence (other than hysteresis) that they 
are important as pointed out by Rygg, et al (l97lt). 

Finally, we consider a simple exanple of the application of this model 
to the interpretation of an observed "hysteresis effect" between the inten- 
sities at two different energies. To illustrate the essential simplicity 
of the model and its basic features we have restricted both energies under 
consideration to the region T > 300 MeV/nucleon where energy loss effects 
are not severe and the simple solution of Parker ( 1963 , 1965 ) remains an 
excellent approximation. Therefore, we can write the stationary solution 
for the density of a given species characterised by magnetic rigidity R and 
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velocity c3 at a given time t simply as 
U(R,3,t) = 


(26) 


tjtnder the assumption that K(R,3,r,^) = K^(r »t )f (R ,6) » where is the 
diffusion coefficient at some reference value of ^ P* Here rj(t) incor- 
porates all the time variations of the modulating parameters and in the 
particular case that is independent of r, we have n(t) = V(t)R(t }/K^(t). 
So long as this assumption of separability holds it follows directly for 
all R and 3 for which equation (26) is valid that 


in 


‘Ul(Rl,3i,t) 

_Ui(Ri,6i,t^) 


Hn 


^2(^2962*'^) 

U 2 (^ 2»&2 


(2T) 


where A = f (R2»62)/^‘(^1»0 i) ^ constant, independent of time. From Equation 27 

we see that under these conditions one expects to be a single valued function of U 
In Figure 5 5 the intensity of protons from 260 to T20 MeV (average 
energy 'v 5OO MeV) observed on a series of balloon flights of the same 
instrument between I965 and 1972 (Rygg and Earl, 1971 » Rygg» 0’ Gallagher 
and Earl, 197^) is plotted versus the Deep River neutron intensity which 
monitors cosmic rays of an average energy 'x. lO GeV. The observations are 
not consistent with the expected single valued relation illustrated by the 
solid curve in the figure which is determined from equation (27) with the 
parameters adjusted so that the curve lies approximately half way between 
the points before solar maximum (solid) and those after (open). The effect 
of a true time lag between the intensities at these two energies has been 
calculated assuming (l) the time variation between maximum and minimum in 
the Deep River intensity is sinusoidal with an eleven-year period and 
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(2) the stationary solution at any time for the 500 MeV proton intensity 
corresponds to the solid ciirve in Figinre 5- Thus the expected proton inten- 
sity for a given neutron intensity including the hysteresis lag given by 
equation (l5) is foxind from the point in the solid curve corresponding to 
the Deep River intensity at a time t earlier under assumption (l). Such 
"loops" are shown superposed on the data in Figure 5 for values of t - 90» 
180, 270 days. The separation between the rising and falling portions of 
the loops can be interpreted as being produced by a value of t between 
180 and 270 days. Of course, this will be so only if there has been no 
time variation in f(R,g) some of which may of course occur (and in fact it 
is just such variations which are proposed to explain hysteresis in all 
stationary models). If we interpret x strictly in terms of this model 
neglecting all other sources of "hysteresis" and higher order effects 
(variations in x over the solar cycle, etc.) discussed elsewhere, it yields 
an estimate of R. With V = 300 Kilometers /second and K(500 MeV protons) 

~ 4 X 10^^ cm^/sec (Jokipil and Coleman, 1968) we find from equation (23) 

R “ V5-55 a.u. which, although somewhat larger than some earlier estimates 
( 0 * Gallagher , I968, Simpson and Wang, 1970) is consistent with current ideas 
and with the small radial cosmic ray gradients observed to date on Pioneer 
10 (Lentz, et al, 1973, Teegarden, et al, 1973, Van Allen 1972). More 
detailed calculations based on this model, incorporating the effects of 
energy loss with computer generated solutions for the complete stationary 
equation have been undertaken and will be reported. The simple example 
above, however, illustrates the utility of the model. 

In c<mtrast to stationary models we see that the model developed here 
provides a physical basis from which a "hysteresis effect" emerges naturally 
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without, the introduction of any new independent modulation parameters. The fact that 
the phase parameter is not independent of E, V and K, thus provides 
a powerful tool for the study of the modulation process itself. Observed 
values of hysteresis phase lags provide direct measures of R and K which 
are not attainable in models based on stationary solutions. 

The author expresses thanks to Drs. T. Rygg and F. Ipavich for several 

helpful discussions. 

This work was supported in part by NASA Grants NGR 21-002-006 and 


NGL 21-002-033. 


26 


APPENDIX A 


We wish to determine the effects of small changes in the modulating 

parameters X, V, K with time by means of a perturbation analysis of the 

modulated density expression in the form of Equation 4, To do this we ex-- 

press the probability distribution function, which is implicitly a slowly 

varying function of time, in terms of its form evaluated at t = t which we 

o 

write as ^ small perturbation 

P(X,V,K,t) - P^(t) + <5P^(t) + A.l 


where 

o 

i 3x At dV At 9K At 

*-0 O O 

in which 


A. 2 


AX AV , AK 
— . — and — 


At ’ At 


are the slow rates of change in each of the modulation parameters. Substi- 
tuting from equations A.l and A. 2 into Equation (4) we obtain 

•'-CO O 


The first term on the right hand side is simply by definition the stationary 

solution U (t=t ) evaluated at t = t . The second term can be evaluated in 
so o 

this first order perturbation treatment by requiring simply that 


P 6t 


<< 


1 


T 


O X 

Physically this Is equivalent to the assumption that most of the contribution 
to the integral takes place on time scales which, although non-zero, remain 
short compared to the scale of the 11-year modulation. Clearly this is a 
better assumption than to assume that 
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1 

P 6t 


= 0 


which is implicit in all stationary solutions. With this assumption the 
derivative with respect to t^ in equation A. 3 can be taken outside the inte- 
gral and the second term becomes : 


^ = XT' h (t )P°U P (t -t)dt| A. 4 

to ' o o o I dtQ I X oJ_„ o o' o I 


from the definition of f ^ in equation 6. Finally this term can be written as 




4t ) 

X o 6 o 


neglecting changes in t (t ) with time since these are second order effects. 

X o 

Sijbstituting back into equation A. 3 one obtains: 


= ^s<^o> k- 


A. 6 


= U (t - T ) 
s o X 

if 5U /6t is approximately constant over the time t . 

Since the exact dependence of P^(t) on the modulated parameters X, V, K 
was not used in any of the above treatment, the result given in equation A. 6 
is quite general and will be valid in the context of any modulation model. 

In particular with a slight change in notation we can perform a perturbation 
analysis on the three dimensional model where the general form for the time 
dependent solution is 

U(R,V,K,r,t^) U^P(R,V,K,r,t^-t)dt A.7 

ico 

analogous to equation 4 in the one dimensional model and P(R,V,K,r,t) is 
given by equation 10. Thus equation (12) and (13) in the text follow 
directly from the above arguments. 
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APPENDIX B 


The hehavior of x is given in terms of integrals of the form 
3c 

r t-n/2 ^-8/t -,t 

•b 

where n is an integer and g and y coefficients depending on X, V, 
and K. Definite integrals of this form diverge in general if either g 
or Y is 0 but otherwise they converge to the solution. 

f t^“^ e dt = 2(|)^ K (2V^) B.l 

•& y y 

(Gradshetyn and Ryzhik, 19^5 » p. 3^0 ) where is a generalized Bessel 

function of imaginary argument. 

Consider first the normalization integral 

_ (yt+x )^ 

/“ P(t)dt = /“ e dt . B.2 

o (UTrKr 'o 

Equation B.l with v = -ht g = X^/Uk, andy= V^/UK yields 


P(t)dt = 


_^/2K 

(UttK)^ 


-VX/2K 

XV ® 


-VX/K 

e 


B.3 
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which can be evaluated in terms of equation B.l with 3 = X^/4K and y * 
(V^/4K + K/X^) and yields equation (l6) in the text directly. 
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APPENDIX C 


The solution for diffusive propagation and escape from a finite 
convecting region in both one and three dimensions can he analyzed in 
terms of the classical probability density w(x,t)(or w(r,t)) following 
the notation of Parker (1965). Consider just a finite one dimensional 
diffusing region with free escape boundaries at x =* ±X and a symmetric 
convective wind of velocity |v| originating at x = o and directed in 
the -X direction for x < o and in the +x direction for x>o. To be 
directly analogous to the physical modulation problem, particles can 
be introduced at both boimdaries so that the entire problem is completely 
symmetric about the origin. 

Thus the diffusion convection equation (2) for w(x,t) is. 


3w(x,t ) _ K 3^w(x,t) V 9w 
is to be solved in the domain o < Ic < 


C.l 

X subject to the boundary conditions 


and 


v(X,t) = 0 
^ (o,t) = 0 


X 

2 

Define the variables s = V t/K and ^ = 
so that equation C.l is separable into 
s’ (s)+ios(s) = 0 
P" U) - P' (?) + (oP(?) = 0 



C.2 


and set w(x,t) = S(s) • L(?) 


C.3 


where u is the separation parameter. The general solutions of C.3 com- 
patible with the boundary conditions are. 


w(x 


,t) = ? Cncos(b ?) • e^^^ e"*^^ 

n=o n 


C.4 


with (On - ^ \ and b = (n + to satisfy C.2. The coefficients 

n n V A 
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Cn are detemiined from the initial condition 

w(x,0) - 6 [x - (X-h)] * C.5 

where h is the distance inside the boundary where the first scattering 
takes place (h « x) by multiplying C.k X(at t = 0) and C.5 by e”^ cos{b^5) 
and integrating from c = 0 to ? = YJ/K yielding 


_ 2 VX/2K ^ 1 \ P < 

= — e cos L(n + ^)ir C.o 

2 VX/2K f , , , vn,, , 

H e ^ (n + ^)7T (-1) /h\ 

Parker (1965) has shown that h can be expressed as h = Uk/v, where v is 
the particle velocity so that the complete solution for the probability 
density P(x,t) = ^ • w(x,t) is given by 


_V(X-x) ^ 

PTX, x,V,K,t ) = ^ e ^ ^(n + h) (-l)’^cos 

X ^ n=0 


(n + Xg)n x 


• exp 


(n + ?s)^Tr%t V ^t 
2 ” ItK 


X 


C.7 


The solution for the similar problem of radial diffusion-convection in 
a three dimensional region inside a free escape boundary of radius E has 
been treated in some detail by Parker (I965) and has the form given 
by equation (12). However the particular radial solutions cannot 
be analytically described as was possible in' C.k above for the one dimen- 
sional case. However the behavior in the limits of large and small K are 
easily treated. 


The solution in the diffusion limit (v»’ 0) is 


P(E,K,r,t) 


2ttK 

Rr 


^ n^ir^Kt 

E (-1) n sin (-^) e ^ 
n=l . ^ 


C.8 
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which is given by Parker ( 1965 ) in a sli^tly different notation. 
Substituting in equation (l3) this yields 

«, n^ir^Kt 

i^CR.K.r) = ^ ^ E 

After evaluation of the integral on the right hand side, equation C.9 
converges to 

Tj.(R,K,V 0) - ^ 

in the limit r « R. C.IO 

Behavior in the limit of small K can be examined in terms of the 
radial diffusion— convection (equation 11) which when written out explicitly 
takes the form 


- If 4 . / 2K ... au V 


C.12 


where the energy loss term has been left but since we are considering 
only spatial propagation. If we further make use of the condition that 

U 9t Tj. 

appropriate to the perturbation analysis in Appendix A the gradient can 

be approximate by its stationary value au vu 

8r " ^ 

so that in the limit K « Vr the terms ^ (3U/3r) and VU/r are small 

r 

compared to V 3U/9r in equation which then takes the approximate 
form 


3U 0. 92u vau „ , , 

— K -r — C.13 

3t 5^2 3r 

which is identical to the one dimensional expression (Equation 2 ) so that 

equations 16 through 19 apply with appropriate substitutions. 

( 
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Figure Captions 

Figure 1. Schematic representation of a one dimensional analog of solar 
modulation. A particle A crosses the boundary (x = X) of a diffusive- 
convective region at t = 0 and is detected by an obsei^rer at x = 0 at a 
later time t. The problem is analyzed in terms of P(XjV^K,t), the probability 
per unit time for observing the particle at x = 0 between given t and t+dt. 
Figure 2. The relative probability distribution as a function of time in 
the "past" (t < t^) for a particle seen now ,(t = t^) at the origin is 
plotted for four sets of "interplanetary" conditions; a) X = Uo a.u. , 

V = 300 Wsec, K = 2 X 10^^ cm^/sec; b)x = UO a.u., V = 300 km/sec, 

K = 4 X 10^^ cm^/sec; c) x = 20 a.u., V = 600 km/sec, K = k x 10^^ cm^/sec, 
and d) X = 20 a.u. , V = 300 km/sec, K = 4 x 10^^ cm^/sec. The gradual 
decay produced by escape from a finite region characterized by a spatial 
scale length X is indicated by the dashed lines. 

Figure 3. The average propagation time delay in a one dimensional finite 
region of scale size X is illvistrated as a function of K and X. The delay is 
proportional to X^/K at large K and is proportional to X/V for small K. The 
/ inlicate values of calculated from the series for the exact solution 
for particles observed at the origin in a symmetric region with boundaries at 
! X. 

Figure k . In a three dimensional radial convective modulating region, the 
average propagation time from a boxmdary R to an observer at r « R can be 
analytically evaluated only in the limits of large and small K indicated by 
the solid lines where the behavior is similar to the one dimensional case. 

The behavior for intermediate K has been represented by a smooth function 
(dashed lines) which joins the two limiting regimes and has a form identical 
to that which is known to be an excellent approximation in one dimension. 
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Typical times of the order of a few months are predicted for some values of 
R,V and K. 

Figure 5 . The "hysteresis effect" observed for 260-720 MeV protons with 
respect to the Deep River neutron intensity. The data is from Rygg and 
Earl, 1971 and Rygg et al, 1974. The behavior before solar maximum (in 
1969) (solid points) and after solar maximum (open points) is not consistent 
with a single valued relation (heavy solid line) but can be interpreted as 
resulting from time lag in the response of the protons of between I8O and 


270 days. 
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